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Abstract. An "exact" method for scalar one-dimensional hyperbolic conservation laws is pre- 
sented. The approach is based on the evolution of shock particles, separated by local similarity 
solutions. The numerical solution is defined everywhere, and is as accurate as the applied ODE 
solver. Furthermore, the method is extended to stiff balance laws. A special correction ap- 
proach yields a method that evolves detonation waves at correct velocities, without resolving 
their internal dynamics. The particle approach is compared to a classical finite volume method 
in terms of numerical accuracy, both for conservation laws and for an application in reaction 
kinetics. 



1. Introduction 

In this paper, a special class of numerical methods for scalar hyperbolic conservation laws 
in one space dimension is presented. An important area in which such problems arise is the 
simulation of nonlinear flows in networks. An example is the flow of vehicular traffic on highways 
|16j . The flow along each network edge is described by a hyperbolic conservation law (e.g. the 
Lighthill-Whitham model [22 for traffic flow). The edges meet at the network vertices, where 
problem specific coupling conditions are imposed (such as the Coclite-Piccoli conditions [18, 3 J for 
traffic flow). Here, we focus on the evolution of the flow along a single edge. In network flows, 
high requirements are imposed on the numerical method. On the one hand, the approach must 
guarantee exact conservation (no cars must be lost), no spurious oscillations must occur (otherwise 
one may encounter negative densities), and shocks (traffic jams) should be located accurately. On 
the other hand, in the simulation of a large network, only very few computational resources can 
be attributed to each edge. 

A commonly used approach is to approximate the governing conservation law by traditional 
finite difference [20] or finite volume methods [13] . Low order methods are generally too diffusive, 
thus do not admit an accurate location of shocks. High order methods, such as finite volume 
methods with limiters [23], or ENO [2J/WEN0 [23] schemes admit more accurate capturing 
of shocks. However, this comes at the expense of locality: stencils reach over multiple cells, 
which poses challenges at the network vertices. Alternative approaches are front tracking methods 
[T7] , These do not operate on a fixed grid, but track shocks explicitly. Thus, shock are located 
accurately. However, smooth parts, such as rarefaction fans, are not represented very well. Another 
class of approaches is based on the underlying method of characteristics. An example is the CIR 
method [3J, which updates information on grid points by tracing characteristic curves. Thus, it is 
a fixed-grid finite difference method. A fully characteristic approach was presented by Bukiet et 
al. PQ that tracks the evolution of particles. Where the solution is smooth, particles follow the 
characteristic curves, and where these curves collide, shocks are evolved. By construction, shocks 
are ideally sharp. While the tracing of the characteristics is high order accurate, the location of 
shocks is only first order accurate. 

Another approach was presented by the authors [HI [1111 [TT] . I n contrast to previous methods, 
here shocks are resolved by the merging of characteristic particles. This is made possible by the 
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definition of a suitable interpolation, which is a similarity solution of the underlying conservation 
law. Hence, in we suggest to call the approach rarefaction tracking. The method admits second 
order accurate location of shocks. In Sect. [2] the fundamentals of this approach, in particular the 
similarity interpolation, are outlined. 

A generalization of the characteristic particle method presented in [10 , namely shock particles, 
is introduced in Sect. [3j A shock particle is a moving discontinuity that carries two function 
values. If the jump height is zero, a classical characteristic particle is recovered. Even though 
shocks are evolved explicitly, the aforementioned similarity interpolation still plays a crucial role, 
thus the approach is fundamentally different from traditional shock/front tracking methods. In 
fact, the presented approach solves the considered hyperbolic conservation law exactly, up to the 
integration error of an ODE. Hence, we call it an exact particle method. In Sect. |4j the evolution 
and interaction of shock particles is shown to give rise to an actual computational method. The 
key idea is that the original partial differential equation is reduced to an ordinary differential 
equation, which then can be solved using a high order ODE solver. A comparison of the exact 
particle approach with a traditional finite volume method (using the package CLAWPACK [2]) is 
presented in Sect. [5] 

The presented approach is highly accurate for hyperbolic conservation laws, and is quite 
amenable to extension to balance laws. Of particular interest here are stiff reaction kinetics, 
in which reactions happen on a faster time scale than the nonlinear advection. In Sect. [6| the 
problem is introduced and some properties of its solution are described. In Sect. [7J a specialized 
adaptation of the particle method is presented. It is based on the exact particle method introduced 
in Sect. [4j with a fundamental adaptation that uses the similarity interpolation to provide a cer- 
tain level of subgrid resolution near the detonation wave. This method is able to track detonation 
waves correctly, without resolving them explicitly Computational results for this application are 
presented in Sect. [8] 



2. Characteristic Particles and Similarity Solution Interpolant 

Consider a scalar conservation law in one space dimension 

«t + (/(«))« = , u(x,0) = u (x) . (1) 

The flux function / is assumed to be twice differentiable and either convex (/" > 0) or concave 
(/" < 0) on the range of function values. We consider an approximation to the true solution of 
by a family of functions, defined as follows. 

Consider a finite number of particles. A particle is a computational node that carries a (vari- 
able) position Xi, and a (variable) function value Mj. Let the set of particles be defined by 
P = {(xi, Ui), . . . , {x n , u n )} with x\ < ■■■ < x n . On the interval [xi,x n ], we define the inter- 
polant Up{x) piecewise on the intervals between neighboring particles, as follows. If u$ = Ui+i, 
then on the interpolant on [xi, ajj+i] is constant Up{x) — Uj. Otherwise the interpolant on [xi, 2^+1] 
satisfies 

x-Xj = f'(Up(x))-f'(ui) 
Xi+i-Xi f(u i+ i) - f'(Ui) 
This defines the inverse interpolant x(Up) explicitly on [xi, 2^+1] . Since / is convex or concave, the 
interpolant Up(x) itself is uniquely defined. As shown in [TU], the interpolation Up defined above 
is an analytical solution of the conservation law ([I]), in the following sense. Consider particles 
moving "sideways" according to P(t) = {(x% + f{u\)t, u{), . . . , (x n + f (u n )t,u n )}. If at time 
t = 0, the particles satisfy X\ < ■ ■ ■ < x n , then for sufficiently short times t > 0, the interpolant 
?7p(t), defined by ^ is the analytical solution to the conservation law at time t, starting 
with initial conditions uo = Up(0). This follows from the fact that each point (x(t),u(t)) on the 
solution moves according to the characteristic equations of ([I]), which are 

( x = f'(u) 
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Figure 1. An illustration of characteristic particles moving according to pj. The 
solution develops a shock at t = T2. The dotted line is the interpolation between the 
particles at time t, the dash-dotted line-after a short time At. The solid line shows the 
solution at t = T2 when the solution develops a shock between particles 2 and 3. 



From the definition of P(t) it is obvious that the particles satisfy ([3]). Any other point is given 
by the above defined interpolation. Replacing xi by xi + f'(ui)t in ([2| and differentiating with 
respect to t yields that x(t) = f'(Up(x)). 

The solution between neighboring particles is a similarity solution that either comes from a 
discontinuity (if the particles depart) or becomes a shock (if the particles approach each other). 
Hence, the solution Up{x) is composed of rarefaction waves and compression waves. Therefore, as 
described in [11], the approach can be interpreted as "rarefaction tracking", which expresses its 
similarity and fundamental difference to front tracking approaches (TTJ [TS] • 

The interpolant Up ^ is a solution of ([!]) until the time of the first collision, i.e. the moment 
when two neighboring particles share the same x-position. For a pair of neighboring particles 
(xi,Ui) and (2^+1, Ui + i), the time of collision is given by 

T = Xl+1 - Xi (A) 

/W) - /'(O ' [> 

If particles depart from each other (i.e. /'(%) < /'(u^+i)), one has T, < 0, thus no collision happens 
in future time. For a set of n particles, the first time of collision is T* = min ({X^ : T t > 0} U 00). 
The solution is continuous until that time. At t = T*, a shock occurs (at Xi, between u^+i and 
i*i), and the method of characteristics alone does not yield a correct solution further in time. An 
illustration of the particle movement and the development of a shock can be seen in Fig. [TJ 

Remark 1. We assume that one is interested in single- valued weak entropy-solutions, which possess 
shocks. In some applications multi-valued solutions are sought, and those can be obtained easily 
by continuing to move the particles according to the characteristic equations Q beyond the 
occurrence of shocks. 

As presented in previous papers [HI OH E] , the use of the method of characteristics even in the 
presence of shocks can be made possible by a suitably designed particle management. Particles 
that collide are immediately merged into a single particle, with a new function value u chosen such 
that the total area under the interpolant Up is preserved. After this merge, the defined interpolant 
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is again continuous, and one can step further in time using solely characteristic particle movement 
([3]). This approach introduces a small error around shocks. Suitable insertion of new particles 
near shocks (before merging) guarantees that the error remains localized near shocks. 



3. Shock Particles 

In this paper, we present an approach that does not introduce any errors intrinsically. While 
the approximation of general initial conditions by finitely many particles involves an error (see 
Sect. 4.1 ), the actual evolution under the conservation law ([I]) is exact — not just pointwise on the 
particles, but in the sense of functions. The new approach generalizes characteristic particles to 
shock particles. 



3.1. Evolution of Shock Particles. A shock particle is a computational node that carries a 
(variable) position Xi, a (variable) left state u~ , and a (variable) right state uf , which satisfy the 
Oleinik entropy condition [7J f'{u~) > f'(uf)). Whenever a shock particle (xi,u~ ,uf) violates 
this conditions (e.g. because it is placed in the initial conditions), it is immediately replaced by two 
particles {xi.uj ,u~) and (x i: uf ,uf), which then depart from each other (since f'{u~) < f'(uf)). 
Given n shock particles P = {(xi, u± , uf), . . . , (x n , u~, uf)}, the interpolant Up on [xi,x r , 
defined piecewise: on [x,-,x,-.fx], it satisfies 

x-Xi f(Up(x))-f(uf) 



IS 



x i+ i - Xi f'(u i+1 ) - f'(uf) 



(5) 



The velocity of a shock particle is given by the Rankine-Hugoniot condition [7] as ii — s(u i , uf), 
where 

S(U ' V) = \ }{u) u = v (6) 
is the difference quotient of /, continuously extended at u — v. 

Remark 2. When implementing this function numerically, one should avoid calculating the dif- 
ference quotient when the distance \u — v\ is very small. For those cases one should consider the 
limiting value f'(u) as a more accurate alternative, or even better, the next order Taylor expansion, 

s(u,u + e)nf'(u) + ±ef"(u), (7) 

can be used. 

At a shock, the function values v7 and change in time as well. Their rate of change is 
exactly such that the interpolation §5§ near the shock evolves as a smooth function should evolve 
under ([T]). Here we derive the evolution for the right value of the shock uf. The argument for 
works analogously. If we had = , the particle would have a function value that is constant 
in time uf = 0, and move with velocity ii = f'(uf). The interpolant ^ with these definitions 
for ii and uf is the correct solution between Xi and Xi + i. If u~ ^ uf , the shock moves at a speed 
different from f'(uf). In order to preserve the same interpolation, the function value uf has to 
evolve according to 

Here x t — f'(uf) is the relative velocity of the shock to a characteristic particle velocity, and 
/ / ("i ) 1 j g ^ s j p e Q £ t j^ e interpolant ([5| at Xi, found by differentiating (|5| with 
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respect to x. The law of motion for a shock particle is thus 

_^ />ti)-/>n i 



( s Kr,u+)-/'Kr)) 



Xi-i-Xi f"{uT) (9) 

/'(T + i)-f«) i 
a*+i - f"(uf) 

Observe that the evolution of a shock particle depends on the neighboring two particles. See Fig. [2] 
for an illustration of the derivation of these equations. 

In the case u~ = uf, we call a shock particle characteristic. In fact, a characteristic particle, as 
described in Sect.[2j is nothing else than a shock particle with jump height zero. This is motivated 
by 

Lemma 1. The motion of a shock particle ([9| reduces to the motion of a characteristic particle 
as uf - ur 0. 

Proof. By definition ([6]), the first equation in (|9| clearly reduces to ±i — f'(ui). In the second and 
third equation, the fraction remains finite while the quantity in the parentheses converges to zero, 
and thus the whole expression vanishes ii~ — >■ 0, and — > 0. □ 

Remark 3. In a numerical implementation, the right hand side (|9| can almost be implemented as it 
stands, with the only modification that the difference quotient (p| is replaced by the characteristic 
speed if \uf — u~ | is less than a sufficiently small value. 

Theorem 2. If the time evolution (JoJ) of the particles P{t) is solved exactly, then for sufficiently 
short times t, the resulting evolution of the interpolation Upu\ is the unique weak entropy solu- 
tion of the conservation law ([!]) with initial conditions uq = U p( ) , on the domain of definition 
[x (t),x n (t)]. 

Proof. Due to the first equation in ([9]), all shocks (including those of height zero) move at their 
correct speeds. By construction, every shock satisfies the entropy condition. Discontinuities 
that violate the entropy condition immediately become rarefaction waves. The second and third 
equation in ^ ensure that each point on the interpolation between particles moves according to 
the characteristic equations □ 

3.2. Interaction of Shock Particles. After some time, neighboring shock particles may collide, 
i.e. share the same x-position. In this situation, the two shocks become a single shock, as the 
following Lemma shows. 

Lemma 3. Two neighboring shock particles (xi,uj~ ,uf) and (xi + i,u~ +1 ,uf +1 ) satisfy — u^ +1 
at their time of collision, if at least one of them is not characteristic. 

Proof. Due to ([9]), the difference in function values between the two shocks evolves according to 
d (..- ,.+\ _ / ^( M r+i."i~+i)-/>r+i) , f(,u+)- s (u7,u+) \ ( f i( - \ f(„.+\\ i 

dt \ U i+l U % ) - I f"(u- +1 ) + /"(«+) J \ J \ 1 W >) x z + 1 -x t 

s(u~ +1 ,uf +1 )-f'(u- +1 ) | f'(uf)-s(u~,u+) \ fii ie \ U i+l ~ U t / 1n \ 

Here £ is a value between u~ +1 — uf given by the Mean Value Theorem. Due to the Oleinik entropy 
condition, both numerators inside the large parentheses are non-positive, and by assumption at 
least one is strictly negative. The signs of /" inside and outside the large parentheses cancel each 
other out. Hence the right hand side always has the opposite sign than u^ +1 — uf. If we assume, 
by negation, that u~ +1 — uf remains finite as Xi+i — Xi — > we get a clear contradiction, and thus 
the difference u~ +1 — uf goes to zero. □ 
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Figure 2. The shock moves with speed s, given by the Rankine-Hugoniot condition, 
which is different than the characteristic speed. Thus shock particles must have varying 
function values. In the figure • denotes moved particles, and Z is the slope of the 
similarity solution to the left of particle 2. 



In the computational method, two shock particles (xi, u i , uf) and (scj+i, u i+11 uf +1 ) that collide 
Xi — Xi + i are simply merged into a single particle (xj, , uf +1 ). Due to Lemma |3j this merge 
does not change the actual solution. Hence, Thm. [2] extends to allow particle merges, assuming 
that the time evolution ^ is integrated exactly. If both interacting particles are characteristic, 
this approach automatically creates a shock. 



4. An "Exact" ODE Based Method 



Due to Thm. [2j the presented approach yields the exact weak entropy solution of the conserva- 
tion law ([I]), when starting with an initial condition uq that can be represented by finitely many 
particles Pq, i.e. uq = Up a . Hence, we call this approach an exact particle method. In practice, two 
types of approximation are performed. First, a general initial function uq cannot be represented 
exactly using finitely many particles, and thus needs to be approximated. This aspect is briefly 
addressed in Sect. 4.1 Second, in general the time evolution ([9| can not be integrated exactly. 



Instead, a numerical ODE solver has to be used. This aspect is addressed in Sect. |4.2 



4.1. Approximation of the Initial Conditions. Whenever the initial function uq can be rep- 
resented exactly by an interpolation Up , one should do so if the number of particles required is 
computationally acceptable. A particular advantage of the presented particle approach is that 
discontinuities can be represented exactly. If the initial function cannot be represented exactly, 
it must be approximated. It is shown in |10| that the interpolation ^ approximates piecewise 
smooth initial conditions with an error of 0(h 2 ), where h is the maximum distance between parti- 
cles, if discontinuities are represented exactly. Furthermore, since the method docs not require an 
equidistant placement of particles, adaptive sampling strategies should be used, such as presented 
in [TU]. These results are based on the particles be placed exactly on the function uq. More general 
approximation strategies that do not have this restriction are the subject of current research. 

4.2. Integration in Time. The characteristic equation ^ can easily be integrated exactly. 
Therefore, characteristic particle movement incurs no integration error, and the next collision 
time between characteristic particles is explicitly given by Q. The particle approach presented 
and analyzed in [9l[T0l[Tl] relies on these properties. The downside of those methods is an intrinsic 
error near shocks. In contrast, the shock-particle method presented in the current paper does not 
incur any errors around shocks. The downside is an error due to the integration of the ordinary 
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Figure 3. Time evolution at t £ {0,0.3,0.6, 1} of the solution to the conservation law 
{I) with f(u) = \u i , both by CLAWPACK and by our particle method. 



differential equation (jjjj. However, it is comparably simple to integrate systems of ODE with very 
high accuracy. In contrast, the construction of high order numerical approaches that approximate 
the PDE directly (such as ENO [Ti]/WENO [23"] or Godunov schemes with limiters [25]), is 
much more challenging. The numerical error analysis shown in Sect. [5] seconds this. 

Another feature (besides accuracy) that the used numerical ODE solver needs to possess is 
event detection. Since at particle collisions, the system undergoes a discontinuous change (the 
number of particles is reduced) , the ODE solver must detect such events with high accuracy. One 
way to do this is to use a solver that can provide a high order interpolation. Several such solvers 
have been derived by Dormand and Prince in [5] for the Runge-Kutta family of ODE solvers. In 
Matlab, event detection is implemented in particular in ode23 .m and ode45 .m. As stated in [23] , 
the latter contains an unpublished^ variation of the interpolation presented in [5] . 

To enhance the performance of the adaptive ODE solver it is helpful to have an estimate 
about the next occurrence of a particle collision. A simple estimate may be obtained by using 
only the first equation of ([9]), thus estimating the collision time between neighboring particles by 

Ti 



s{u 



Xj+i- 
:+i» u i+i 



) — s(u i ,u^) 



5. Numerical Error Analysis of the Particle Method 



We investigate the order of accuracy of the presented particle method, and compare it to the 
benchmark PDE solver CLAWPACK [2J [2T] . The comparison between a finite volume method 
and a particle method is tricky, since the two approaches are fundamentally different. First, the 
finite volume approach works with average function values in fixed cells, while with particles the 
interpolation {T5J) defines a function everywhere. This difference can be overcome by constructing 
errors in the Lr sense from cell averages, as described in [lOj . Second, the finite volume method 
has a fixed spacial resolution Ax, while particles move, merge, and are generally anything but 
equidistant. Third, in a convergence analysis of a finite volume method, the spacial resolution 
and the time step are chosen proportional At = C Ax. In contrast, the particle method becomes 
exact if At — > 0, assuming that the initial conditions can be represented exactly by finitely many 
particles. 

Here, we consider an initial condition that can be represented exactly by the interpolation 
The reason is that we do not want to measure the error in approximating general initial conditions 
(for this aspect, please consult PU]). Instead, we want to investigate the error in the particle 
evolution. We consider a second order and a fourth order accurate Runge-Kutta method for the 
time evolution of ([9]). Times of particle collisions are found and resolved with the same order of 
accuracy. For the CLAWPACK runs, we specify a desired CFL number [4] of 0.8 and let the code 
choose At as it finds suitable. In practice, for this problem, this amounts to having At s» Ax. 

Specifically, we consider the conservation law ( Jip with flux function f(u) — \u A : and initial 
function uq(x) = Up (x), which is the interpolation (121) defined by the characteristic particles Pq = 
{(0, 0.1), (0.1, 0.1), (0.2, 0.9), (0.4, 0.9), (0.5, 0.7), (0.6, 0.7), (0.7, 0.1), (1.0, 0.1)}. The time evolution 



Details can be found in the Matlab file ntrp45.m. 
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of the solution is shown in Fig. [3] in four snapshots at t g {0, 0.3, 0.6, 1}. In fact, what is shown 
is the solution obtained by the particle method, integrated with an accuracy that the error is not 
noticeably in the eye norm. For a comparison, we show the results obtained by a second order 
CLAWPACK method, with Ax = 0.05. 

The convergence of the error for various approaches is shown in Fig. [4] We consider the 
i 1 ([0, 1]) error at time t = 1. One can observe that the overall order of the particle method 
equals the order of the ODE solver used. Thus, with the standard RK4 solver, machine accuracy 
is obtained already for moderate time steps. Note again that we do not need to increase the 
number of particles to obtained convergence. This is a crucial advantage of the presented particle 
method. In comparison, the CLAWPACK method yields an order of convergence slightly better 
than first order. This is particularly due to the presence of the shocks in the solution, and the 
large derivatives in the initial condition. 



6. Stiff Reaction Kinetics 



Many problems in chemical reaction kinetics can be described by advection-reaction equations, 
for which the reaction happens on a much faster time scale than the advection. We consider the 
balance law 

t*t + (/(«)) x =#0, (11) 

where < u < 1 represents the density of some chemical quantity. The advection is given by the 
nonlinear flux function /, which is assumed convex (the case of / concave is analogous) and to 
be of order 0(1). The reactions are described by the source tp. Here, we consider a stiff bistable 
reaction term 

^(«) = i«(l-u)(u- i 9) (12) 
where < j3 < 1 is a fixed constant. This source term drives the values of u < (3 towards 0, 
and values of u > /3 towards 1. The reactions happen on a much faster time scale 0(t), where 
t«1. This example is presented for instance in [15] . Since the source term tp does not act in a 
discontinuity, equation (111 possesses shock solutions as the homogeneous problem ([!]) does. In 
addition, it has traveling wave solutions that connect a left state ul ~ with a right state ujjwl 
by a continuous function. To find those solutions, we make a traveling wave ansatz u(x, t) = 
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where £ = is the self similar variable, 
differential equation for v 

= 



This transforms equation (111 into a first order ordinary 



For v = (3, the numerator of (131 vanishes. A solution that connects a state Vl < ft to a state 



vr > (3 can only pass through v = if the denominator of (131 vanishes as well, i.e. r = /'(/?), 
which yields the velocity of the traveling wave. The shape of the wave is then found by 

-rt> 



integrating (13) using r = /'(/?)■ Since u(x,t) — v(^P L ), the traveling wave has a thickness 
(in the ^-coordinate) of 0(t). This analysis of traveling wave solutions is in spirit similar to 
detonation waves of reacting gas dynamics |12j . The value /3 plays the role of a sonic point in 
detonation waves. The advection-reaction equation (11) is studied in [8]. The traveling wave 



solution ( 13 ) results from a balance of the advection term (which flattens the profile) and the 
reaction term (which sharpens the profile). Since r is very small, these traveling waves look 
very similar to shocks, yet they face the opposite direction, and travel at a different velocity (if 

/'(/?) * /(i) - /(o)). 

In computations, the recovery of the exact shape of the traveling waves is typically not very 
important. However, the recovery of their correct propagation velocity is crucial. As described 
in |21j . equation (11) can be treated in a straightforward fashion using classical finite volume 



approaches. However, correct propagation velocities of traveling waves are only obtained if these 
are numerically resolved. Thus, with equidistant grids, one is forced to use a very fine grid 
resolution h = 0(t), which is unnecessarily costly away from the traveling wave. This problem can 
be circumvented using adaptive mesh refinement techniques, however, at the expense of simplicity. 
An alternative approach, presented in [15] . yields correct traveling wave velocities even with grid 
resolutions h 3> 0(t), by encoding specific information about the structure of the reaction term 
into a Riemann solver. 



7. A Particle Method for Stiff Reaction Kinetics 

Here, we present an approach based on the particle method introduced in Sect.|4]that uses the 
"subgrid" information provided by the interpolation ^ to yield correct propagation velocities of 
traveling waves, without specifically resolving them. The characteristic equations for (|11[) are 



(14) 



As before, our goal is to generalize these characteristic equations to obtain an evolution for shock 
particles. This requires the definition of an interpolation. We use the interpolant as if there 
were no reaction term. Clearly, this is an approximation, and the resulting method is not exact any- 
more. At any time t, we define the solution by shock particles P(t) — {(xi, u^), . . . , (x n , u~, m+)}, 
and the interpolation C/p(t), defined by (J5|. Adding the reaction term to Q, we now let the par- 
ticles move according to 




/ K )) 



i 



Xi-1 - Xi /"(it. ) 
m)-/'K + ) 1 



f'(u 



/"«) 



) 



Hut) 



(15) 



where the shock speed s(u i , uf) is defined as before by ([6]). By construction, shocks move at their 
correct velocity, and for a characteristic particle 1*7" = uf , ( 15 1 reduces to the correct characteristic 



evolution (14). Clearly, this approach does not remove the stiffness in time. Hence, an implicit 



ODE solver should be used. System (15) yields an accurate solution on the particles themselves, 



as well as an accurate evolution of shocks. However, traveling waves, as given by (13) are not 



represented well. The reason is that each particle moves very quickly towards or 1. Then, the 
reaction term is not considered anymore, since ip(0) = and ^(1) = 0. In order to correctly 
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Figure 5. Correction approach for the advection-reaction equation with dominant re- 
action term. The vertical dashed lines denote the three roots of the source term. 



represent traveling waves, the continuous solution that goes through the sonic point ft has to be 
considered. We do so by the following correction approach. 

We assume that the bistable nature of the reaction term and the value of the unstable root 
ft are known. Whenever the solution increases (in x) from a value u < ft to a value u > ft, a 
special characteristic particle is placed at u = ft, that moves with velocity x — f'(ft)- We call 
such a particle a sonic particle. Each particle that neighbors a sonic particle is treated in a special 
way. As a motivation, consider the situation shown in Fig. [5] A left state is followed by a sonic 
particle, which is followed by a right state 1. The interpolant shown is ^ for f(u) — \u 2 . The 
thick solid graph shows the initial configuration at time t. The dotted graph shows the obtained 
solution when evolving the particles for a time step At according to the method of characteristics 



(14 1. Since -0 vanishes on all particles, the source is neglected. Hence, this approach does not lead 



to correct traveling wave solutions. The thin solid graph shows the correct evolution of the initial 
function, considering the interpolation. This function cannot be represented exactly by particles 
and the interpolant ^ , but it can be approximated by modifying the x- values of the two particles 
that neighbor the sonic particle, in such a way that the areas under the solution both left and 
right of the sonic particle are reproduced correctly. The dashed graph shows the function that 
results from this correction. Below, we describe the correction approach in detail. 

7.1. Computational Approach. Consider a particle u~_ 1 ,uf_ 1 ) that is a left neighbor of 

a sonic particle {xi,ft). The case of a right neighbor particle is analogous. Let the interpolant on 
[xi-i,Xi] be denoted by U(x), and its inverse function X{u). From §5§ it follows that X'(u) = 
. f"(u). Using the interpolant U(x), we can integrate the reaction term between the 

two particles. The substitution rule yields 

ip(U(x))dx= / ijj(u)X'(u)du = 







ij)(u)f" (u)du . 



(16) 



np)-nutx)Juu 

This expression represents the full influence of the reaction term on the continuous solution between 
the value uf_ 1 and the sonic value ft. As derived in [10] . the area under the interpolant on [xi-\, Xi] 
is given by 



U (x)dx = (xi — Xi-i) a(u 



i-l'l 



where a(v,w) = — ^p^p is a nonlinear average. Now consider a new particle (xi-\ 
Ai£i_i, uf_ 1 , u~l_i) be inserted between x^\ and Xi. This insertion changes the area by 

f(ft){ft - ut x ) - (f(ft) - /(tijjti)) 



AA = Axi-i (a(w+ 



= A 



Xi-l- 



nft)-nuU) 



(17) 
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Figure 6. Time evolution at t £ {0.1,0.2,0.3,0.4} of the advection reaction equation 
( 11 I with t = 0.01. The thick gray graph shows the true solution, while the dots denote 



the particle approximation. 



Equating the rate of area change given by (17), and expression (16), yields 

where c(v, w) = f , {w )(J^)-utw)-f(v)) ■ The scaling / = 0(1) and ip = 0(±) implies that c(v, w) = 
O(-), if to — v = 0(1). A similar derivation for the right neighbor yields that a new particle 
(xi + i — Axj+i, u^ +1 , u^ +1 ) needs to be inserted with 



At 



c ( u i+l,P) ( x i - x i+l) ■ 

Due to the bistable nature of the reaction term, one will frequently encounter a nearly constant 
left state, i.e. both \uf_ 1 — uJ_A and \u^_ 2 — u^_A are very small. In this case, the particle i — 1 
can just be moved by Axi-i, instead of creating a new particle. Using a characteristic particle 
notation only, the resulting modified evolution equations are 

±i-i = f'(ui-i) + c(ui_i,/9) (xi - Xi-i) 

Ui-l = ^(Ui-l) 
Xi = f(P) 
Ui = 

x i+i = f'(lk+i) + c(u i+1 ,/3) (x l+ i 
Ui+i = ip(v,i+i) ■ 



< 



Xi) 



This implies that 



^ (x 



dt v -« = Cf'O 9 ) - /'(^i-l)) - cK-1,/3) (Xi -Xi-i) , 

i.e. the distance Xi — x,_i converges to an equilibrium value - ■ Since c(uj_i, /3) = OQp), 

the equilibrium distance between the sonic particle and its neighbors is 0(t). Hence, the presented 
approach yields a traveling wave solution, represented by three particles that move at the correct 
velocity /'(/?), and whose distance from each other scales, correctly, like O(r). 



8. Numerical Results on Reaction Kinetics 

For assessing our method numerically, we compare it to the benchmark PDE solver CLAW- 
PACK [2], for the advection-reaction equation ( [Tl"] ). We consider the reaction term (12) with 
(3 = 0.8, and choose the Burgers' flux f(u) = ^ir. Four different values for the reaction time 
scale are considered: r e {0.1,0.024,0.008,0.004}. The spacial resolution is Ax = 0.02. As initial 
condition, we use u(x, 0) = 0.9 exp (— 150(a; — \) A \ For solving this problem using CLAWPACK 
we simply use the code from the CLAWPACK website [3TJ Chapter 17]. This code was written 
specifically to solve this stiff Burgers' problem. 

The time evolution of the solution of ( 11 1 is shown in Fig.[6]in four snapshots at t € {0.1, 0.2, 0.3, 0.4}. 
The thick grey graph shows the true solution. The solid dots denote the particle method. One can 
see that at t = 0.1 the solution is still in the transient phase, since characteristic particles are still 
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Figure 7. Computational results for the advection reaction equation (111 with r 6 
{0.1,0.024,0.008,0.004}. The thick grey graph shows the true solution, while the dots 
denote the particle approximation. 



visible on the (soon-to-be) detonation wave. At t = 0.2, the wave structure is almost converged. 
The plots at t = 0.3 and t = 0.4 show how the detonation wave catches up to the shock. 

Figure[7]shows the solution at the final time t = 0.4, for four choices of r £ {0.1, 0.024, 0.008, 0.004}. 
The thick grey graph shows the true solution. The solid dots denote the particle method. The 
circles show the CLAWPACK results. Note that for r = 0.1, the solution is still in the transient 
phase, while for the other values, the detonation wave is comparably well established. For the 
selected resolution Ax = 0.02, CLAWPACK successfully captures the shock for r e {0.1,0.024}. 
The detonation wave for t = 0.024 is nicely represented as well. However, CLAWPACK clearly 
fails to resolve the shock and the detonation for r — 0.004. The intermediate value r = 0.008 is 
on the edge of failure. In comparison, the particle method works for all values of r. The shock 
is optimally sharp, and the detonation wave moves at the correct velocity and has the correct 
width. The trouble that CLAWPACK is having with these equations is due to the stiff source. 
The problem is that the width of the shock is always O(Ax), but this is too large when r becomes 
small. The source is too active both in the detonation shock and in the regular forward-facing 
shock. This leads to incorrect shock speeds. 



9. Conclusions and Outlook 



We have presented a particle method that solves scalar one-dimensional hyperbolic conservation 
laws exactly, up to the accuracy of an ODE solver and up to errors in the approximation of the 
initial conditions. The numerical solution is defined everywhere. It is composed of local similarity 
solutions, separated by shocks. A numerical convergence analysis verified this accuracy claim 
for the flux function f(u) = u 4 /4. In this example, the basic RK4 method yields solutions up 
to machine accuracy using a few hundred time steps. Since general initial conditions can be 
approximated with second-order accuracy (see [TU]), the overall method is at least second order 
accurate, even in the presence of shocks. 

The method has also been extended to balance laws that describe stiff reaction kinetics. The 
tracking of a sonic particle in combination with a correction approach for neighboring particles 
yields a method that evolves detonation waves at correct velocities, without actually resolving 
their internal dynamics. The evolution of the sonic particle comes naturally in the considered 
particle method, while for classical fixed grid methods, a similar approach is much less natural. 
Numerical tests show that the particle method approximates the true solutions very well, even for 
fairly stiff systems, for which CLAWPACK fails due to an under-resolution of the wave and the 
shock. 

The philosophy of the considered application in stiff reaction kinetics is that one can find 
efficient approaches for more complex problems by using the exact conservation law solver at the 
basis. It is the subject of current and future research to apply the same philosophy in other 
applications. Examples are nonlinear flows on networks. The presented particle method can be 
used to solve the actual evolution on each edge exactly. While an approximation has to be done 
at the network nodes, it is plausible to believe that this approach yields more accurate results 
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that classical method that are far from exact on the edges themselves. Further generalizations to 
consider are the treatment of higher space dimensions using dimensional splitting, and systems of 
conservation/balance laws. 
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